{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fracture extrusion\n",
    "\n",
    "This tutorial shows how to create a 3D fracture network from a 2D outcrop. The idea is to provide a more complete documentation, an example, and highlight assumptions and simplifications in the module."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What can this module do?\n",
    "1. Create a 3D fracture network that is consistent with a flat 2D outcrop, in the sense that if the 3D network is cut at the outcrop height, the 2D network should be reproduced.\n",
    "2. Produce circular fracture planes, potentially cut, with a prescribed or random dip angle.\n",
    "3. Preserve Y/T-intersections, in the sense that if one fracture ends in another in the outcrop, this relation will be kept in the extruded network. However, fractures that meet outside the outcrop plane will not be cut in this way.\n",
    "\n",
    "\n",
    "## What is missing \n",
    "1. A geological foundation for the way the extrusion is carried out. What is here so far is an attempt to apply common sense, but nothing more.\n",
    "2. Population of fractures that do not intersect with the outcrop plane. These will have to be assigned e.g. by stochatsic rules, but this is beyond the scope of PorePy.\n",
    "3. Extrusion to other than circular shapes.\n",
    "4. ...\n",
    "\n",
    "### Will any of the missing pieces be implemented shortly?\n",
    "There are no concrete plans for this.\n",
    "\n",
    "## What is the intended use for the model?\n",
    "Simulations in quasi-realistic (somewhat generously interpreted) geometries. Typically, a verification that developed numerical schemes etc. can handle first test towards real geometries.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How to use\n",
    "\n",
    "First, import statements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import porepy as pp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The starting point for the extrusion is an outcrop, described by endpoints and connection between points.\n",
    "\n",
    "Note that fractures should *not* be split into branches before applying extrusion. The reasons are \n",
    "1. The extrusion is carried out by drawing circles through the endpoints of the fractures (the outcrop trace is a chord). Lines through branches will create all sort of problems.\n",
    "2. Abutting relations (T-intersections) are identified by comparing original and split (by the extrusion algorithm) traces."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAADO9JREFUeJzt3G+IZXd9x/HPZ3bHiZddEWqgks3c\na0GkaaCbZggrgc6wBLpGsVQoKKOPhHnSQgRBDPug9cGlzyRP+uTGhBbmoggKLakiETMjQic6Gyey\ncQxE2RkXpYsUcZeV6SR+++Be9+5m5+78OSfzm/me9wsu7D177jm/+eXe956cc+Y6IgQAyGOi9AAA\nAPUi7ACQDGEHgGQIOwAkQ9gBIBnCDgDJ1BZ22yds/9j2C3VtEwCwf3UesT8lab3G7QEADqCWsNs+\nI+mjkr5Sx/YAAAd3sqbtPCPpC5JOj1vB9oKkBUn6E+nRzqlTNe0aABrg97/Xpbfe+k1E3L/bqpXD\nbvtjkq5FxCXbc+PWi4iepJ4kzZw+HavXr1fdNQA0x9ycvLy8sZdV6zgV87ikj9u+Iulrks7bXqxh\nuwCAA6gc9oh4OiLORERH0iclfS8iPl15ZACAA+E+dgBIpq6Lp5KkiFiStFTnNgEA+8MROwAkQ9gB\nIBnCDgDJEHYASIawA0AyhB0AkiHsAJAMYQeAZAg7ACRD2AEgGcIOAMkQdgBIhrADQDKEHQCSIewA\nkAxhB4BkCDsAJEPYASAZwg4AyRB2AEiGsANAMoQdAJIh7ACQDGEHgGQIOwAkQ9gBIBnCDgDJEHYA\nSIawA0AyhB0AkiHsAJAMYQeAZAg7ACRD2AEgGcIOAMlUDrvt+2z/0Partl+z/aU6Bobm6Pf76nQ6\nmpiYUKfTUb/fLz0k4Fg7WcM2tiSdj4gbticl/cD2tyNipYZtI7l+v6+FhQXdvHlTkrSxsaGFhQVJ\n0vz8fMmhAcdW5SP2GLgxfDo5fETV7aIZLl68eCvqf3Tz5k1dvHix0IiA46+Wc+y2T9hek3RN0osR\n8fIO6yzYXrW9ur29XcdukcDm5ua+lgPYXS1hj4i3IuKspDOSHrP98A7r9CJiJiJmJicn69gtEpie\nnt7XcgC7q/WumIj4raQlSRfq3C7y6na7arVadyxrtVrqdruFRgQcf3XcFXO/7fcO//xuSU9I+lnV\n7aIZ5ufn1ev1NDU1JUlqt9vq9XpcOAUqqOOumPdL+nfbJzT4h+LrEfFCDdtFQ8zPz+vZZ5+VJC0t\nLZUdDJBA5bBHxE8kPVLDWAAANeA3TwEgGcIOAMkQdgBIhrADQDKEHQCSIewAkAxhB4BkCDsAJEPY\nASAZwg4AyRB2AEiGsANAMoQdAJIh7ACQDGEHgGQIOwAkQ9gBIBnCDgDJEHYASIawA0AyhB0AkiHs\nAJAMYQeAZAg7ACRD2AEgGcIOAMkQdgBIhrADQDKEHQCSIewAkAxhB4BkCDsAJEPYASAZwg4AyVQO\nu+0Hbb9ke932a7afqmNgAICDOVnDNt6U9PmIeMX2aUmXbL8YET+tYdsAgH2qfMQeEb+OiFeGf74u\naV3SA1W3CwA4mFrPsdvuSHpE0ss7/N2C7VXbq9vb23XuFgBwm9rCbvuUpG9I+lxE/O7tfx8RvYiY\niYiZycnJunYLAHibWsJue1KDqPcj4pt1bBMAcDB13BVjSc9JWo+IL1cfEgCgijqO2B+X9BlJ522v\nDR9P1rBdAMABVL7dMSJ+IMk1jAUAUAN+8xQAkiHsAJAMYQeAZAg7ACRD2AEgGcIOAMkQdgBIhrAD\nQDKEHQCSIewAkAxhB4BkCDsAJEPYASAZwg4AyRB2AEiGsANAMoQdAJIh7ACQDGEHgGQIOwAkQ9gB\nIBnCDgDJEHYASIawA0AyhB0AkiHsAJAMYQeAZAg7ACRD2AEgGcIOAMkQdgBIhrADQDKEHQCSIewA\nkAxhB4Bkagm77edtX7N9uY7tNUG/31en09HExIQ6nY76/X7pIeEI4H0xwlxUEBGVH5L+WtJfSbq8\nl/UfPXUqmmxxcTFarVZIuvVotVqxuLhYemjFzM7OxuzsbOlhFMX7YoS52MHsbEhajT001jEIc2W2\nO5JeiIiHd1t35vTpWL1+vZb9HkedTkcbGxt3LZ+amtK5c+cKjKi8tbU1SdLZs2cLj6SclZUVbW1t\n3bW83W7rypUrhz+ggsZ9Rpo4F7fMzcnLy5ciYma3VQ/tHLvtBdurtle3t7cPa7dH0ubm5o7Ld/pQ\noznG/fcf937JbNzP3MS5OIiTh7WjiOhJ6kmDI/bD2u9RND09PfZoZGlp6fAHdATMzc1JUmN/fmn8\nUer09HSB0ZQ17jPSxLk4CO6KKaDb7arVat2xrNVqqdvtFhoRjgLeFyPMRTWEvYD5+Xn1ej21223Z\nVrvdVq/X0/z8fOmhoSDeFyN/nIupqSlJavRcHEQtF09tf1XSnKT3SfofSf8UEc+NW7/pF09xN07F\nYCe8L26zj4untZxjj4hP1bEdAEB1nIoBgGQIOwAkQ9gBIBnCDgDJEHYASIawA0AyhB0AkiHsAJAM\nYQeAZAg7ACRD2AEgGcIOAMkQdgBIhrADQDKEHQCSIewAkAxhB4BkCDsAJEPYASAZwg4AyRB2AEiG\nsANAMoQdAJIh7ACQDGEHgGQIOwAkQ9gBIBnCDgDJEHYASIawA0AyhB0AkiHsAJAMYQeAZAg7ACRT\nS9htX7D9uu03bH+xjm1m1+/31el0NDExoU6no36/X3pIxfT7fa2srGh5eZm54H1xC++LCiKi0kPS\nCUk/l/Rnkt4l6VVJD93rNY+eOhVNtri4GK1WKyTderRarVhcXCw9tEPHXIwwFyPMxQ5mZ0PSauyh\ny45BnA/M9ocl/XNE/M3w+dPDfzD+ZdxrZk6fjtXr1yvt9zjrdDra2Ni4a/nU1JTOnTtXYETlrKys\naGtr667lzMVIu93WlStXDn9ABY37jDRxLm6Zm5OXly9FxMxuq9ZxKuYBSb+87fnV4bI72F6wvWp7\ndXt7u4bdHl+bm5s7Lt/pQ53duJ+ZuRgZ937JbNzP3MS5OIiTNWzDOyy7638DIqInqScNjthr2O+x\nNT09PfZoZGlp6fAHVNC9jsyYi4Hp6ekCoylr3GekiXNxEHUcsV+V9OBtz89I+lUN202r2+2q1Wrd\nsazVaqnb7RYaUTnMxQhzMcJcVLSXE/H3emhw1P8LSR/Q6OLpX9zrNU2/eBoxuDjUbrfDdrTb7UZf\nFGIuRpiLEebibQ7z4qkk2X5S0jMa3CHzfETc85/Vpl88BYB928fF0zrOsSsiviXpW3VsCwBQDb95\nCgDJEHYASIawA0AyhB0AkiHsAJAMYQeAZAg7ACRD2AEgGcIOAMkQdgBIhrADQDKEHQCSIewAkAxh\nB4BkCDsAJEPYASAZwg4AyRB2AEiGsANAMoQdAJIh7ACQDGEHgGQIOwAkQ9gBIBnCDgDJEHYASIaw\nA0AyhB0AkiHsAJAMYQeAZAg7ACRD2AEgGcIOAMkQdgBIhrADQDKVwm77722/ZvsPtmfqGhQA4OCq\nHrFflvQJSd+vYSwAgBqcrPLiiFiXJNv7e+GNG9LcXJVdA0CzrK3tedVKYd8P2wuSFoZPt7y8fPmw\n9n3EvU/Sb0oP4ohgLkaYixHmYuRDe1lp17Db/q6kP93hry5GxH/sdTQR0ZPUG25zNSI4Jy/m4nbM\nxQhzMcJcjNhe3ct6u4Y9Ip6oPhwAwGHhdkcASKbq7Y5/Z/uqpA9L+i/b39njS3tV9psMczHCXIww\nFyPMxcie5sIR8U4PBABwiDgVAwDJEHYASKZY2Jv+dQS2L9h+3fYbtr9Yejwl2X7e9jXbjf7dBtsP\n2n7J9vrws/FU6TGVYvs+2z+0/epwLr5Uekyl2T5h+8e2X9ht3ZJH7I39OgLbJyT9q6SPSHpI0qds\nP1R2VEX9m6QLpQdxBLwp6fMR8eeSzkn6hwa/L7YknY+Iv5R0VtIF2+cKj6m0pySt72XFYmGPiPWI\neL3U/gt7TNIbEfGLiPg/SV+T9LeFx1RMRHxf0v+WHkdpEfHriHhl+OfrGnyIHyg7qjJi4Mbw6eTw\n0dg7PWyfkfRRSV/Zy/qcYy/jAUm/vO35VTX0A4yd2e5IekTSy2VHUs7w1MOapGuSXoyIxs6FpGck\nfUHSH/ay8jsadtvftX15h0djj06HdvrWtMYejeBOtk9J+oakz0XE70qPp5SIeCsizko6I+kx2w+X\nHlMJtj8m6VpEXNrra97RLwHj6wjGuirpwduen5H0q0JjwRFie1KDqPcj4pulx3MURMRvbS9pcB2m\niRfYH5f0cdtPSrpP0ntsL0bEp8e9gFMxZfxI0gdtf8D2uyR9UtJ/Fh4TCvPg+6+fk7QeEV8uPZ6S\nbN9v+73DP79b0hOSflZ2VGVExNMRcSYiOhq04nv3irpU9nbHg34dwbEXEW9K+kdJ39HgAtnXI+K1\nsqMqx/ZXJf23pA/Zvmr7s6XHVMjjkj4j6bztteHjydKDKuT9kl6y/RMNDoRejIhdb/PDAF8pAADJ\ncCoGAJIh7ACQDGEHgGQIOwAkQ9gBIBnCDgDJEHYASOb/ASHRRQ5PncqOAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f8f604bd9e8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Define points and connections\n",
    "pt = np.array([[0, 0], [2, 0], [1, 0], [1, 3], [0, 1], [2, 1], [3, 0], [3, 1]]).T\n",
    "edges = np.array([[0, 1], [2, 3], [4, 5], [6, 7]]).T\n",
    "\n",
    "# Temporary min-max coordinates, only valid for 2d domain\n",
    "domain = {'xmin': -1, 'xmax': 4, 'ymin':-1, 'ymax': 4}\n",
    "pp.plot_fractures(domain, pt, edges)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The edges are defined so that fractures 0 and 1 define a Y-intersection, 1 and 2 an X-intersection, while no 3 is isolated.\n",
    "\n",
    "Next, we extrude the fractures, collect them in a FractureNetwork, and export to Paraview."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fractures = pp.extrusion.fractures_from_outcrop(pt, edges)\n",
    "network = pp.FractureNetwork(fractures)\n",
    "network.to_vtk('extrusion_1.vtu')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The resulting network looks as follows\n",
    "\n",
    "<img src=\"fig/extrusion_1.png\" alt=\"Drawing\" style=\"width: 400px;\"/>\n",
    "\n",
    "\n",
    "Observations:\n",
    "1. The Y-intersection between 0 and 1 is preserved in the 3d version of the network, while fractures 1 and 2 cuts through each other.\n",
    "2. Although fracture 3 was isolated in the outcrop, it cuts fracture 0 somewhere underneath.\n",
    "\n",
    "Also note that, as the height of the fracture center was decided randomly by the extrusion, the network will look different from one realization to the next.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Assumptions\n",
    "1. All extruded fractures will be circles, for lack of more information.\n",
    "2. The fractures are assumed to be cut at an arbitrary height relative to their center. In the extrusion, the z-coordinate of the center, thus the radius of the fracturce, is drawn as random numbers. The behavior can be overruled by the parameter 'extrusion_angle', which should be between 0 and $\\pi$. Note that the limiting cases define infinite fractures, while $\\pi/2$ will give a point contact in the outcrop plane, none of which are physical.\n",
    "3. Depending on the fracture radius, the above procedure may yield a T-intersection where the terminated fracture extends beyond the constraining fracture. To avoid this, the size of the constraining fracture will be increased, that is, the z-coordinate of the center is modified. This behavior can be overruled by setting the optional parameter 'ensure_realistic_cuts' to False.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# T-intersections and tolerances\n",
    "Whether two lines are found to intersect in a T, X, or not at all depends on the tolerance used in geometric calculations. While this applies to most geometry computations in PorePy, there is a particular issue that may arise when working with outcrop data: When fractures are collected into a database, points in T-intersections are routinely snapped to the abutting fracture. In the decision to snap, a tolerance is applied, at least implicitly.\n",
    "\n",
    "When the extrusion model is invoked, it will identify T-intersections based on a user defined tolerance that may not be the same as in the snapping process. Moreover, it is not clear that the two tolerances mean the same. We have experienced some issues with what from the user side was considered T-intersections was treated as X-intersections by the extrusion module because of this.\n",
    "\n",
    "To deal with this:\n",
    "1. Inspect the extruded fracture network, and see if you find the expected T-intersections.\n",
    "2. If necessary, redo the snapping with a tolerance sufficient to make extrusion behave well. The function snap_points_to_segments in the comp_geom model can be used for this.\n",
    "\n",
    "It should also be quite easy to rewrite the extrusion model to accept a pre-computed topology, however, this has not been prioritized."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# More advanced features\n",
    "In addition to the basic features described above, the extrusion module has some more advanced features that goes further towards creating more realistic networks. These are however also less tested. It is not clear whether the approach will ever create something resembling geological realism. Anyhow, the features are:\n",
    "1. Assign inclines to angles, so that they are not all vertical.\n",
    "2. Use family-based parameters, so that fractures that are grouped together (say, were created under the same geological regime) are assigned similar parameters, but different from other families."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
